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ABSTRACT 

We investigated the instability of advective accretion flow as a consequence of angular momen- 
tum transfer in one-dimensional, quasi-spherical transonic accretion flow around a non-rotating 
black hole. The code is designed to include the effects of viscosity; the hydrodynamics com- 
ponent preserves angular momentum strictly with Lagrangian and remap method in absence of 
viscosity, while the viscosity component updates viscous angular momentum transfer through the 
implicit method. We performed two tests to demonstrate the suitability of the code for accretion 
study. First, we simulated the inviscid, low angular momentum, transonic accretion flow with 
shocks around a black hole, and then the subsonic, self-similar ADAF solution around a New- 
tonian object. Both simulations fitted the corresponding analytical curves extremely well. We 
then simulated a rotating, viscous, transonic fluid with shocks. We showed that for low viscosity 
parameter, stable shocks at larger distance are possible. For higher viscosity parameter, more 
efficient angular momentum transfer in the post-shock disk makes the shock structure oscillatory. 
Moreover, as the shock drifts to larger distances, a secondary inner shock develops. We showed 
that the inner shock is the direct consequence of expansion of the outer shock, as well as creation 
of regions with dl/dr < due to more efficient angular momentum transfer near the inner sonic 
point. We showed that all disk parameters, including emissivity, oscillate with the same period 
as that of the shock oscillation. Our simulation may have implication for low frequency QPOs, 
e.g., GRO Jf655-40 and XTE J1550-564. 

Subject headings: accretion - hydrodynamics - instabilities - methods: numerical 



1. INTRODUCTION 

Investigation of the fiow behavior of the ac- 
creting matter in the vicinity of a black hole is 
important since the spectrum and the intensity of 
the emitted radiation depend on the flow struc- 
ture. The event horizon presents the unique inner 
boundary condition in which the in-falling mat- 
ter crosses the horizon with the speed of light (c). 
Therefore black hole accretion has to be transonic, 
as a result of which existence of one sonic point or 
critical point is assured for black hole accretion. 
General relativity also ensures that matter must 
posses sub-Keplerian angular momentum closer to 



the horizon. Although within the marginally sta- 
ble circular orbit (vms), the angular momentum at 
f ^ fms is deflnitely sub-Keplerian (and the value 
I'^lms — l\r=r^J, but at larger radius the angular 
momentum should be generally large. Therefore, 
a general accretion disk should have viscosity to 
remove the angular momentum outwards. The 
first serious model of viscous ac c retion disk was 
presented bv lShakura fc Sunvaev ( 1973 ). in which 
the angular momentum distribution was Keple- 
rian, the accretion disk was geometrically thin 
and optically thick. In Shakura-Sunyaev disk the 
pressure and advection terms were not properly 
considered and no attempt was made to satisfy 
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the inner boundary condition around a black hole 
apart from the adhoc termination of the disk at 
T < Tms- Along with this theoretical short com- 
ing, Shakura-Sunyaev disk also failed to explain 
power-law high energy part of a black hole candi- 
date spectrum. Therefore, search for another com- 
ponent of an accretion disk which may explain the 
origin of the high energy radiations from black hole 
candidates, were undertaken by various groups. 
One such mo del which g o t a w id e attention was 
ADAF [e.g., Ilchimaru l (|l977[) . iNaravan fc Yi 



(|1994| ) hereafter NY94]. This model was first 
constructed around a Newtonian gravitational po- 
tential, where the viscously dissipated energy is 
advected along with the mass, momentum and 
the entropy of the flow. The original ADAF solu- 
tion was self-similar and wholly subsonic, and was 
found to be thermally and dynamically stable. 
Howover, the low viscosity ADAF showed con- 
vecti ve instability ( Igumenshchev fc Abramowicj 
19991) . that has no dynamical effect if the an- 



gular momentum is transported outward but it 
is dynamically important in case the opposite is 
true. The global solution of ADAF showed that 
the flow actually becomes transonic at around few 
Schwarzschild radii (r^), and the self-similarity 
may be maintaine d far away from the sonic point 



( Chen Rt al\\im% 



Simultaneous to these developments, there were 
some interesting research going on sub-Keplerian 
flows around black holes. It has been shown 
that sub-Keplerian flow does posses multiple sonic 
point in a significant range o f the energy-angular 
mom entum parameter space ( Liang &: ThompsonI 
1980l ). One of the consequences of existence of 
multiple sonic points, is that the flow accreting 
through the outer sonic point can be slowed down 
by the centrifugal barrier. This slowed down mat- 
ter acts as barrier to the faster fluid following it. 
If the strength of the barrier is s trong enough 
then accretion shocks may form (jChakrabarti 
19891) . General global solutions in the advec- 



tive domain incorporating viscosity and thermal 
effects were obtained by many independent re- 
searchers (Chakrabarti 1990, 1996; Lu et al. 199^ 
iLanzafame et aLlll998l : IGu fc Lu 1120041 ). Further- 
more, it has also been shown that the global 
ADAF solutio n is a subset of the general advec- 

Il999l ). Whether a flow 



tive solutions (Lu et al. 



brid solution with or without shock will depend 
on the outer boundary condition and the physical 
processes dominant in the disk. 

Although steady-state solutions are possible in 
a cer t ain range of parameter spac e dChakrabarti 



1989: Chakrabarti k Das I 12004 iMolteni et al. 



1994 Il996a[ ). but advective solutions with dis- 



continuities such as shocks are generally prone 
to various kind of instabilities. Since, various 
flow variables across the shock surface jumps 
abruptly, this results in a markedly different cool- 
ing, heating and other dissipation rates across 
the shock. This may render the shock unsta- 
ble. For example, in presence of bremsstrahlung 
cooling, resonance between cooling timescales and 
in fall timescales in the post shoc k part of the 



disk gi ves rise to oscillating shoc ks (jMolteni et al. 



Il996bh . .Lanzaf ame et al. ( 1998h showed that be- 
yond a critical viscosity post-shock disk may os- 
cillate. The interaction of the outflow and the in- 
flow may also cause the b en ding instability in the 



disk ( Molteni et a/.M2001l ). IMolteni et al 



]y m tne 
] (|l999l) 



showed that in presence of non-axisymmetric az- 
imuthal perturbations the shock initially becomes 
unstable but stabilizes within a finite radial ex- 
tent into an asymmetric closed pattern. More- 
over, the post-shock region may be associated 
with the elusive Compton cloud that produ c es the 
hard photons ( Chakrabarti fc Titarchuk ' 'l995|: 



ICh a krabarti fc Manda l 2006: Mandal fc Chakarabarti 
20 0^ and may also be the base of the jet 



(Das fc Chakrabarti 'l999':'Das et a;."2001:' Chattopadhvav fc Das" 
5007; Das fc Cb attopadhvav 2008; Becker et all 
[2008t bas et aLl (20091 ). Therefore instabilities of 
the post-shock region may manifest itself as the 
variabilities observed in the emitted hard pho- 
tons se en in microquasars and active-galactic 
nuclei 



aclei Jl 
dst, [Fi 



Molteni et al\ Il996b ). To add a new 



twist, mrkumura fc Tsuruta 



(|2004 ) c onjectured 



the p resence of multiple shocks and iGiri et al. 



will follow an ADAF solution or some kind of hy- 



(l2010[ ) actually reported the presence of two os- 
cillating shocks giving rise to two quasi-periodic 
oscillations. 

In this paper, we concentrate on the study of 
instabilities of rotating fluid around black holes, 
generated by the angular-momentum transport by 
viscosity. Since the temperature, density etc are 
higher and the velocity is lower in the post-shock 
region compared to the pre-shock region, the angu- 
lar momentum transport rate should be different 
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in the two regions of the disk. In other words, in 
this paper we simulate transonic, viscous, rotating 
fluid around black holes. We employ a new code 
to study the effect of angular momentum trans- 
port in the accretion disk. Unlike other purely 
Eulerian codes, this new code is especially devel- 
oped to strictly conserve angular momentum in 
absence of viscosity. In §2, governing equations 
and assumptions are presented. In §3, the code 
which was built to calculate the evolution of an- 
gular momentum as accurately as possible is de- 
scribed, along with tests for a rotating transonic 
flow and a viscous flow. In §4, the structure and 
the instability shown in simulations are presented, 
along with descriptions on the nature of the insta- 
bility. A summary and discussion is presented in 
§5. ' 

2. BASIC EQUATIONS 

The one-dimensional time-dependent equations 
for quasi-spherical accretion of viscous flows are 
given by 



0, 



dVr 

'dt 
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dr 
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1 dp 
p dr 
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$i and e are the gas density, ra- 
dial velocity, specific angular momentum, gravita- 
tional potential and specific internal energy, re- 
spectively. The angular velocity is defined as 
= Z/r^. The suffix i in equation (2) denotes 
N or PN, correspo nding to Newtonian or pse udo- 
Newtonian gravity ( Paczvnski fc Wiita Ill980l) . re- 
spectively, and are given by 



dr 

where p, Vr 
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and 



$PN = -- 



GM, 



BH 



(5) 



(6) 



where Mbh is the black hole mass and the 
Schwarzschild radius is rg = 2GMbh / ■ The 



pseudo-Newtonian potential is widely used to 
mimic the Schwarzschild geometry. For the gas 
pressure the equation of state for ideal gas is as- 
sumed, i.e., 

P^il- l)pe, (7) 

where 7 is the ratio of specific heats. For viscosity, 
the a prescription (Shakura & Sunyaev 1973) can 
be assumed, i.e., the dynamical viscosity coeffi- 
cient is described by 



p ~ ap 



n 



K 



(8) 



where 



TP 
P 



(9) 



is the square of the adiabatic sound speed, and 



O - 



r dr 



1/2 



(10) 



is the Keplerian angular velocity, and the viscosity 
parameter a is a constant which is generally less 
than 1. It is to be noted that the actual expression 
oi^lK depends on the gravitational potential used. 
Finally following NY94, the parameter / measures 
the fraction of the viscously generated energy that 
is stored as entropy and advected along with flows. 
The value f — 1 corresponds to the limit of full 
advection and has been used in this paper. 

In the following, we use c and rg as the units of 
velocity and length, respectively, unless otherwise 
stated. In the geometrical units, the unit of time 
is = rg/c. 

3. CODE AND TESTS 

One of the most demanding tasks in carry- 
ing out numerical simulations of equations (1) - 
(4) is to calculate the evolution of the angular 
momentum as accurately as possible. Capturing 
shocks sharply should also be important in resolv- 
ing structures with clarity, if shocks are involved. 
It has been known that the latter can be achieved 
by using codes based on modern, upwind finite- 
difference schemes on an Eulerian grid. However, 
without viscosity in such Eulerian codes, it is nor- 
mally the azimuthal momentum (pv^) and not the 
angular momentum that is treated as a conserved 
quantity. On the other hand, codes based on the 
Lagrangian concept, such as the SPH code, can 
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be designed to preserve the angular momentum 
strictly. Although it has been successfully applied 
to many studies of accretion flows however, the 
SPH code is known to be unduly dissipative (see, 



e.g. 



Molteni et al. 1996al for discussions) 



Here, we describe an Eulerian code that was 
built to accurately calculate the evolution of the 
angular momentum including its transport due to 
viscosity, and at the same time to capture dis- 
continuities (shocks and contact discontinuities) 
sharply with minimum numerical dissipation. The 
code is composed of two parts; hydrodynamic and 
viscosity part. The hydrodynamic part is based on 
the Lagrangian TVD plus remap approach. The 
Lagrangian/remap approach is not new in numer- 
ical hydrodynamics and was employed previously 



(jColella fc Woodward 1 119841 ). But here we show 



that in this approach the equation for angular mo- 
mentum conservation can be directly solved, so 
the hydrodynamics part can be designed to pre- 
serve the angular momentum strictly in the ab- 
sence o f viscosity. At the same time , the TVD 
scheme ( Harten il983 : Rvu et allll993[ ) guarantees 
sharp reproductions of discontinuities and mini- 
mum numerical dissipation. In the viscosity part, 
the viscous angular momentum transfer is updated 
through an implicit method, assuring it is free 
from numerical instabilities related to it. But the 
viscous heating is updated with a second order ex- 
plicit method, since it is less subject to numerical 
instabilities. 

3.1. Hydrodynamic Part 

The hydrodynamic part consists of the La- 
grangian step and the remap step. First in the 
Lagrangian step, the equations for Lagrangian hy- 
drodynamics are solved. On the Lagrangian grid 
defined with mass coordinate, equations (1) - (4), 
except for the centrifugal force, gravity, and vis- 
cosity terms which are treated separately [see be- 
low], can be written in a conservative form as 



dm 
dp 



0, 



dt 

dt am 



dE d{r^Vrp) 



dt 



dm 



(11) 
(12) 
(13) 
(14) 



where r and E are the specific volume and the 
specific total energy, respectively, that are related 
to the quantities used in equations (1) - (4) as 



1 



E = e 



(15) 



The mass coordinate is related to the spatial co- 
ordinate via 

dm = p{r)r^dr, (16) 



and its position can be followed with 



dr 
'dt 



Vr{m,t). 



(17) 



Equations (11), (12), (14) form a hyperbolic 
system of conservation equations, and upwind 
schemes can be applied to build codes that ad- 
vance the Lagrangian step using the Harten's TVD 
scheme, which is an explicit, second-order, finite 
difference scheme to sol ye a hyperboli c system of 



conse rvation equations (jHarten 1 119831 : iRvu et al. 
19931) . We note that the angular momentum in 
Eq. (13) is preserved, so it need not be updated 
in the Lagrangian step. 

In the remap step, the quantities evolved in the 
Lagrangian grid are redistributed to the Eulerian 
grid, to preserve the spatially fixed grid structure. 
Before the Lagrangian step, the Lagrangian and 
Eulerian grid zones coincide. But after the step, 
the Lagrangian grid zone moves to the updated 
position 

where v is the time-averaged velocity, and so it 
does not coincide with the Eulerian grid zone any 
more. Not only are the quantities conserved in the 
Eulerian grid, the density, radial momentum, and 
total energy, but also the angular momentum are 
remapped. For the remap, we employ the third 
order accurate scheme used in the PPM code (see 
Colella & Woodward 1984, for details). 

With the Lagrangian and remap steps, equa- 
tions (1) - (4) are updated in the Eulerian grid, 
except for the centrifugal force, gravity, and vis- 
cosity terms on the right hand side. The centrifu- 
gal force and gravity terms are calculated sepa- 
rately after the Lagrangian and remap steps such 
that 



^hydro ^lag-frcmap 



+ At 



remap 



dr 



(19) 
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Then the viscosity terms are calculated, as dis- 
cussed in the next subsection. 

Non-uniform Eulerian grids can be employed in 
the code. For the problem in this paper we use a 
grid, where the size of cells increases exponentially 
as 

Ar, = Arix5*-\ (20) 

to achieve higher resolution at the origin. Here 
Ari the size of the first grid cell and S is the in- 
crement factor. 

3.2. Viscosity Part 

Viscosity has two effects on accretion flows. 
First, it transfers the angular momentum out- 
wards, allowing the matter to accrete inwards. At 
the same time, it acts as friction, which results in 
viscous heating. 

Since the term for the angular momentum 
transfer in Eq. (3) is linear in Z, it can be solved 
implicitly. Substituting of (r™ + r<='"'^P)/2 for I, 
Eq. (3) without the advection term becomes 



The evolution of an inviscid flow, which en- 
ters the outer boundary with a small amount 
of angular momentum and approaches a black 
hole described by Paczyi iski & Wiita potential 



( Paczvnski fc Wiita 1980f) . is calculated in cylin- 



i+l 



jrcmap 



rt\ Ticmap jrcmap 



(21) 

forming a tridiagonal matrix. Here a^, hi, and 
Ci are given with p, /x, and r as well as Ar and 
At. The tridiagonal matrix can be solved easily 
for Z"*"^ w ith an appropriate boundary condition 
([Press et al. 1992). The term for the viscous heat- 
ing in Eq. (4) is also linear in e (note that cx e), 
so it alone can be solved implicitly too. However, 
combining the two linear equations for I and e 
becomes complicated. Through numerical exper- 
iments, we found that the explicit treatment for 
the viscous heating does not cause any numerical 
problem. So instead of implementing a compli- 
cated scheme to solve simultaneously I and e im- 
plicitly, we solve the angular momentum transfer 
implicitly, while solving the viscous heating explic- 
itly. 

3.3. Tests 

Two tests are presented to demonstrate that the 
code can handle transonic flow as well as viscous 
flow which are involved in our problem. The ca- 
pability of the code to capture shocks sharply and 
resolve structures clearly is tested with a transonic 
accretion flow. In these tests we reproduce well- 
known results. 



drical geometry. We note that previous subsec- 
tions describe the code only in spherical geome- 
try, but the code is actually written in arbitrary 
geometries. For this test, the version in the cylin- 
drical geometry is used. Without viscosity, the 
angular momentum is preserved. A shock can 
form, if the rotating flow through the outer sonic 
point approaches the centrifugal barrier and de- 
celerates discontinuously, due to the twin effect 
of centrifugal force and pressure. In the test, 
I = l.Scrg, a value slightly below the marginally 
stable value [Ims — (3/2)'^/^crg], is used. The val- 
ues of the flow quantities in geometrical units are 
{p,p,vr) = (0.71809,0.007604,-0.083566) at the 
injection radius ri„j = bOvg. The sound speed at 
finj is Cs = 0.1188c, hence the fluid is subsoni- 
cally injected. The fluid becomes supersonic af- 
ter crossing the outer sonic point at rco = 27.9rg, 
becomes subsonic at the shock r^h — 7.89rg and 
enters the blackhole supersonically after crossing 
the inner sonic point at rd = 2.563rg. For the 
ratio of specific heats, 7 = 4/3 is used. This is 
the case of a s table s tanding shock considered in 
Molteni et all (|l996a[ ). In Figure 1, the numeri- 
cal solutions (open circles) of density p (top) and 
radial Mach number Mr = Vr/cg (bottom) using 
2048 uniform grid cells is compared with the an- 
alytical solution (solid line). The fiow in smooth 
regions coincides with the analytical solution very 
well and the shock position matches the analyti- 
cal value very well. The agreement of analytical 
result with the current code is better than those 
with the purely E ulerian TVD code and the SPH 
code presented in lMolteni et all (1996a). 

Next, the performance of the code for a sub- 
sonic viscous flow is tested with a self-similar 
ADAF solution. Matter is steadily injected with 
Keplerian angular velocity into the computational 
domain at rinj, and the simulation lasts until the 
steady state is reached. Figure 2 presents the flow 
quantities after the steady state is reached and 
compares them with the analytic solution. The 
Newtonian potential is used, and the ADAF is de- 
scribed with a self-similar solution (NY94). The 
values of physical parameters used in this test 
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are 7 — 4/3, and a = 0.3. The simulation was 
performed on an exponentially increasing grid of 
780 cells with An = 0.4972rs and b = 1.01. 
Here is the sink size. The injection position is 
Tinj ^ 3.6915 X lO'^rg. The quantities are drawn in 
units of the Keplerian velocity and the Keplerian 
angular momentum at the sink, and Ixirs), 

and the density is in an arbitrary unit. The fig- 
ure shows that the analytic solution is reproduced 
very closely in the region between r ^ lOr^ and 
r ~ lO^Tg in a box of size 1.1611 x lO^rg. The er- 
ror in the specific angular momentum is less than 
a few % at most. 

4. RESULTS OF SIMULATIONS 

In previous numerical works, oscillation phe- 
nomena in accretion flows around black holes re- 
lated to the QPO were reported. The study of 
inviscid supersonic accretion flows around a New- 
tonian central object, showed that the accretion 
disk with shock struc ture to be dynamically unsta- 
ble (|Rvu et aLlll995h . Global transonic accretion 
flows around black holes have been known to ex- 
hibit stationary shocks for inviscid (*Chakrabarti 



B-989: Molteni et Q/...1996a:lDas a/.i.2001.) as well 
as di s sipative flows ( Chakrabarti Il996t Lu et al. 



I999I: iLanzafame et a/.lll998HBecker et aLll2008l ). 
However, since the post-shock flow is hotter, 
denser and slower, the dissipation rate in the post- 
shock flow is shorter than that of the pre-shock 
flow, which may make the post-shock flow unsta- 
ble. Indeed it has been shown that the energy- 
angular momentum parameter space for stand- 
ing shock d ecreases with the increase of viscosity 
paranieter dChakrabarti fc Das _ 2004 ;_ Gu fc Lu 



2004 iDas et aLll2009l ). ILanzafame et al\ (|l998l ) 
simulated viscous transonic flow and showed 
steady shocks exist for low viscosity, while for 
highe r viscosity shock becomes unstable. How- 
ever, ILanzafame et al\ (|l998) restricted their in- 
vestigations for a hot flow {Tinj ^ lO^^X at the 
injection radius), and very low viscosity parame- 
ter [a ^ fewx 10"'^). Since the pre-shock flow was 
chosen to be hot (post-shock disk was obviously 
even hotter), the angular momentum removal was 
very efficient in both the pre-shock as well as post- 
shock disk, even when the viscosity parameter was 
low. The length scale of the computation box was 
only about few tens of Vg . In this paper we simu- 
late viscous transonic flow which is cold to begin 



with, and investigate the instability arising from 
reasonably higher viscosity of the flow. The rea- 
son to choose cold flow at the injection is to have 
a very different angular momentum transport rate 
in the post-shock and pre-shock disk, and thereby 
to maximize the effect of shock instability. More- 
over, we keep the length scale of computation box 
fairly large so as to study large amplitude and low 
frequency shock instability. 

4.1. Shock Formation in Inviscid Rotating 
Flow 

We start our set of simulations with a low 
energy, rotating, transonic, inviscid flow around 
a black hole (described by ^pn)- The steady- 
state, inviscid, transonic solution corresponds to 
a flow characterized by Bernoulli parameter or 
speciflc energy [which in our unit system is £ ~ 
0.5z;2 + c2/(7 - 1) + ;2/(2r2) - l/{2(r - 1)}] and 
specific angular momentum {I). Parameters £, I 
for the inviscid flow are 1.25xl0~^c^ and 1.8crg, 
respectively. The ratio of specific heat is given 
by 7 = 1.4. The Bondi radius (the length scale 
within which gravity becomes important) is de- 



fined asrs = GMsH/cj^, where i 



f(7-l) 



for an inviscid flow. In this particular case vb = 
10*^rg and Cs,oo = 7.071 xlO-*c. The analyti- 
cal, steady-state solution of flows for these pa- 
rameters gives two physical sonic points, the in- 
ner one is at rd = 2.394rg and the outer one at 
rco = 199991.04rg « 0.2rB. The analytical so- 
lution also predicts a shock at Vgh = 22.2rg. It 
has been shown in connection to Figure 1, that 
it is possible to simulate transonic flow quite ac- 
curately with subsonic injection i.e., when rinj > 
rco- However in the present scenario, rgh ^ rco- 
Therefore, if rinj > fco, then a large amount of 
computation time will be wasted in simulating un- 
interesting region of the disk. Hence, without any 
loss of generality, we choose the injection param- 
eters from the supersonic portion of the analyt- 
ical curve in order to reduce computation time. 
To further reduce the computation time and also 
to achieve higher resolution close to the center, 
we use exponentially-increasing grids, which are 
3,553 cells with Ari = 0.0296 and 5 = 1.001 and 
the length of the computation box corresponds to 
lOOOrg — O.OOlrs. The injection radius is hence 



— lOOOr^ — O.OOlrs, and the flow radial ve- 
locity [vr{ini) — 2.970 x 10~^c], specific angu- 
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lar momentum [linj = l.Scrg] and sound speed 
[cs{inj) = 4.827 x 10"'^ c] at rinj is taken from the 
analytical solution. In Figure 3, we compare the 
steady-state analytical solution (solid) with sim- 
ulation result (open circle) when steady state is 
reached. Various flow variables like p (a), Vr (b), 
Cs (c) and p (d) are plotted with log{r). Figure 3, 
shows excellent agreement of the simulation result 
with the analytical curve and the shock has been 
captured within 2-3 cells. 

4.2. Shock Oscillation of Viscous Flow 

Time dependent solutions of viscous transonic 
accretion flow are obtained by starting with the 
inviscid flow described in §4.1 as the initial condi- 
tion, and then increasing the viscosity parameter 
a. The action of the viscosity can be understood 
from equations (3,4). For accretion {i.e., when 
Vr < 0), if the r.h.s of equation (3) is negative 
then the angular momentum will be transported 
outwards. In presence of Shakura-Sunyaev type 
viscosity, the angular momentum transport in the 
post-shock subsonic flow is much more efficient 
than the pre-shock supersonic flow. As a result of 
which, angular momentum starts to pile up in the 
immediate post-shock fluid, resulting in a jump 
in the angular momentum distribution across the 
shock. Similarly, equation (4) tells us that the vis- 
cous heat dissipation in the post shock disk will 
also be higher compared to the pre-shock disk. 

It is well known that a standing shock forms if 
the total pressure ( ram-|-gas pressure) i s conserved 



across the shock ( Chakrabarti 19891 ). In Fig- 



ure 4a, the inviscid solution produces a stationary 
shock location Vgh — 22.2rg, as is shown in Figure 
3. In Figure 4b, the shock location as a function of 
t is plotted for a = 0.003. The excess gas pressure 
due to viscous heat dissipation, and the increased 
centrifugal force due to the piled-up angular mo- 
mentum in the post-shock disk, pushes the shock 
front outwards. For low a, shock front moves to 
a larger location {vsh ~ 31rg as in Fig. 4b) where 
the balance between total outward push and the 
total inward pressure from the pre-shock flow is 
restored. However, for higher viscosity parameter 
a = 0.006 (Fig. 4c), the enhanced angular mo- 
mentum transport creates an even more stronger 
outward push and the shock front overshoots a 
possible equilibrium position and the shock starts 
to oscillate. Interestingly, when the shock moves 



to around ~ 70 rg and beyond, a second shock 
tends to emerge, which expand and collide with 
the outer shock. The combined shock then drifts 
outwards, the inner shock re-emerges and the cycle 
continues. In the following, let us make a detailed 
investigation on transonic flow with higher viscos- 
ity, and the emergence of two shocks. 

In Figure 5a - 5d, we have plotted radial ve- 
locity Vr (dashed-dot) and the sound speed Cs 
(sohd) at four time steps (a) t — 2.9xl0^rg/c, 
(b) t = 3.165xl05rg/c, (c) t = S.AxlO'^rg/c and 
(d) t = 3.615 xlO^rg/c, where the outer boundary 
conditions are same as Figure 3, and a = 0.01. 
In Figure 6a - 6d, the speciflc angular momen- 
tum distribution l(r) is plotted for the same sim- 
ulation and for the same time steps as in Figure 
5a - 5d. Corresponding panels of Figures 5 and 
6 are to be considered in tandem to understand 
this complicated phenomenon. The different snap 
shots in Figures 5a - 6d, corresponds to (i) the 
maxima in outer shock (Figs. 5a and 6a), (ii) the 
expansion stage of the combined shock just after 
the minima in the next cycle (Figs. 5b and 6b), 
(iii) just before the maxima of the outer shock 
(Figs. 5c and 6c), and (iv) just after the max- 
ima of the outer shock (Figs. 5d and 6d). In 
Figure 5a, there are two shock structures, the in- 
ner shock is at rsh(in) ^ 130rg and the outer 
shock rs;i(out) ~ 500rg. Corresponding angu- 
lar momentum distribution in Figure 6a shows 
that the dl/dr > in the range r < 20rg and 
rs;i(in) < r < rs/i(out), while dl/dr < is in 
the range 20rg < r < rs;i(in). In Figure 5b, the 
two shocks merge and the combined shock is at 
rs;i(in) = Vshiont) = Vsh ~ lOOrg. Figure 6b 
shows that dl/dr > in a region where r < 20rg, 
and dl/dr < for 20rg < r ^ rgh, with a smaller 
hump in angular momentum distribution around 
rsh- The angular momentum distribution attains 
a tall peak, and the enhanced centrifugal pressure 
almost stalls the infall (Fig. 5b) in that region. 
However, due to the extra pressure from the piled- 
up I the combined shock moves outwards, while 
the contact discontinuity wave resulting from the 
collision of shocks moves towards the black hole. 
As a result the sound speed {i.e., temperature) in 
the immediate post-shock ffow drops (Fig. 5c), 
and the angular momentum distribution becomes 
dl/dr > (Fig. 6c). This allows for a more free 
infall and u,. in the immediate post shock disk in- 
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creases considerably (Fig. 5c) . As the contact dis- 
continuity wave is absorbed by the black hole the 
angular momentum in the immediate post-shock 
region gets reduced considerably, becomes su- 
personic in the region. However, the flow again hit 
the centrifugal barrier closer to the black hole and 
the inner shock rcemerges (Figs. 5d and 6d). Wc 
note that the regions oi dl/dr < are subject to 
the rotational instability. The non-steady behav- 
ior shown here should be partly attributed to the 
instability. 

4-2.1. On emergence of the inner shock, shock 
collision and the angular momentum 
transfer 

In the top panel of Figure 7, the shock oscilla- 
tion is plotted for a = 0.01. Therefore, each of the 
panels of Figures. 5a - 6d corresponds to the vari- 
ous snap shots of flow variables taken from the top 
panel of Figure 7 (time sequences a d arc marked 
on the figure). Similar to Figure 4c, rsh in the top 
panel of Figure 7 starts to oscillate as the viscosity 
is turned on. A transient inner shock, i.e., rshijo) 
develops when rs/i(out) ^ 80rg. Initially the dy- 
namics of the two shocks are similar to that of Fig- 
ure 4c, i.e., the inner shock forms when rs/i(out) 
is at the maxima, and then rs/j(in) collides with 
the contracting outer shock, and the merged shock 
then expands. However, for t > 0.2xlO^Tg the 
shock dynamics slowly change, both shock ex- 
pands and then contracts, and the shocks collide 
while contracting. The merged shock then reaches 
a minima and then expands, and this cycle con- 
tinues. The query about the formation of the in- 
ner shock can be understood as follows: as the 
original shock expands to a distance ^ &Org, the 
sound speed in the immediate postshock region 
and close to the black hole differs by almost an 
order of magnitude. Hence the rate of angular 
momentum transport in a region nearer to the 
inner sonic point is much higher than the region 
closer to the shock. Hence, the angular momen- 
tum transport rate is not only markedly different 
between the post-shock and pre-shock region, but 
also within the post-shock flow when the shock 
expands to a very large distance. Hence the an- 
gular momentum piles up in between inner sonic 
point (rci) and the shock {e.g., Fig. 6b), which 
enhances the centrifugal barrier and impedes the 
accretion. Continued shock expansion reduces the 



post shock sound speed (i.e., temperature) and 
creates a mild but positive angular momentum 
gradient, which increases the infall velocity in the 
immediate post-shock flow. This can continue up 
to the extent that the post shock fluid once again 
become supersonic in the immediate post-shock do- 
main, however further down-stream the piled-up 
angular momentum virtually stops the supersonic 
inflow causing the formation of the inner shock. 
The inner shock again jacks up the temperature, 
which causes the inner shock to expand too. If 
the outer shock is contracting then the two shocks 
may collide, or, both the shock may expand in 
phase and collide during the contraction phase. 
The combined shock then expands and the whole 
cycle is repeated. It may be noted that the inner 
shock emerges half way into each of the cycles and 
hence it is a persistent feature. 

It is interesting to seek the radiative property of 
such oscillatory dynamics of the disk. We estimate 
the bremsstrahlung loss aposteriori from the disk, 
as a representative of the radiative loss. It is well 
known that the bremsstrahlung emission cx Cg. 
In Figure 7b, we plot Erotai Br/pfnj Cs{inj) as a 
function of time, where 

Erotai Br= Cs r^dr, (22) 

Pinj and Cs{inj) are the density and sound speed 
of the flow at the outer boundary or rj„j . Hence, 
the bottom panel of Figure 7 represents total 
bremsstrahlung emission from the computational 
box compared to the bremsstrahlung emission at 
rinj. Interestingly, bremsstrahlung emission also 
has a periodic behavior, whose period is similar 
to the period of the shock oscillation. However, 
the shock maxima/minima may or may not co- 
incide with either emission maxima or minima. 
In this particular case, initially there is no dis- 
tinct correlation, but as the oscillation approaches 
quasi saturation, the emission maxima coincides 
with the combined shock minima. And the emis- 
sion minima coincides with the rising phase of the 
combined shock and when the inner shock has not 
been formed. As the combined shock contracts it 
pushes the post-shock matter inward just like a 
'bellow' in a blacksmiths shop. This jacks up the; 
p, Cs and Vr. The enhanced p and Cg contribute 
to form the emission max;ima. As the combined 
shock expands, the flow variables like p, Cg in the 



8 



immediate post region decreases, and the emis- 
sion starts to decrease until it reaches the minima. 
Since the wide difference of Cs also triggers the dif- 
ferential / transfer in the inner regions and outer 
regions of the post-shock disk, the angular momen- 
tum again starts to pile up and starts the forma- 
tion of the inner shock described above. There ex- 
ists a secondary peak in the bremsstrahlung emis- 
sion as well which appears to be related to the 
dynamics of the inner shock. The time lag be- 
tween the shock maxima and the emission maxima 
is 6t ~ 2xlO*Tg. Initially the oscillation period of 
the shock was T^^^ ^ 5 x lO^Tg. The oscillation 
period gradually increases to a quasi-saturation 
value of Tosc 



■ SxlOVg. Since, 



2GM 



10" 



H 



Mr, 



sec. 



therefore. 



Tosc ~ 8 X I0\a ~ 8 X 10-1 



Mr. 



sec. 



(23) 



(24) 



This would correspond to the frequency ~ 
0.125^2; for stellar mass black hole i.e., Mbh ~ 
lOM©. In case of a super-massive black hole 
{Mbh ~ 10* M0), these time-scales will corre- 
spond to 2.5 yr variabilities. However, since there 
are two shocks, we are interested to see the influ- 
ence of the dynamics of the two shocks on emis- 
sion. In the left panel of Figure 8, we plot the outer 
shock (top), inner shock (middle) and the relative 
bremsstrahlung emission (bottom) for reference, 
and in the right panels we plot the power den- 
sity spectra for the outer shock (top), inner shock 
(middle) and the bremsstrahlung emission (bot- 
tom) for a stellar mass black hole [Mbh = IOMq). 
The power density spectrum of the outer shock 
shows a frequency of 0.125Hz. The power den- 
sity spectrum of the inner shock has a prominent 
peak at the frequency ~ 0.25Hz and a weaker 
peak around ~ 0.125i?2;, the secondary peak sug- 
gests that the oscillation of the outer shock forces 
a weak periodicity on the inner shock as well. The 
power density spectrum of the inner shock is a bit 
noisy, since the time variation of the inner shock, 
although persistent, is not continuous. Interest- 
ingly, Figure 8, shows that the bremsstrahlung 
emission also peaks at around the same frequen- 
cies as that of the two shocks, confirming that 



the quasi periodicity in the emission is due to the 
quasi periodicity of the two shocks. 

In presence of such a dynamical disk, it is in- 
triguing to investigate the time variation of the 
amount of matter and angular momentum con- 
sumed by the black hole. Let us define the mass 
loss parameter, or the ratio of the rate of mass 
cannibalized by the black hole to the rate of mass 
injected, as M/Minj, where M = {pvrr'^)\sink and 
Minj = {pVrr^)\rinj ■ The angular momentum loss 
rate is defined as L/Linj, where L = Ml sink and 
Minjlinj- The average specific angular mo- 



mentum of the disk is defined as 

fldr 
fdr 



< I >= 



(25) 



In Figure 9, we plot the mass-loss parameter 
(top panel), angular momentum loss rate (mid- 
dle), and the average angular momentum of the 
disk (bottom) as a function of time. The pro- 
files of the mass-loss parameter and the angular 
momentum loss rate are similar to that of the 
bremsstrahlung emission rate. Since the distribu- 
tion of p peaks when the shock is at the minima, 
the peaks of the mass-loss parameter and the an- 
gular momentum loss rate coincide with the peak 
of the bremsstrahlung emission. As the shock re- 
cedes, Vr and p decreases resulting in matter get- 
ting accumulated in the disk i.e., M/Minj < 1- 
As the shock contracts it squeezes more matter 
into the black hole (accumulated in the expan- 
sion phase) than it is being supplied, therefore 
M/Minj > 1- It is to be noted that in this particu- 
lar case, the disk prefers to stay in the state where 
M/Minj ^ 1- The angular momentum loss rate 
follows M / Minj ■ Interestingly, the maxima of the 
average angular momentum of the disk coincides 
with the minima of the emission, mass-loss rate 
and the angular momentum loss rate. The bot- 
tom panel of Figure 9 suggests that if the average 
angular momentum of the disk increases, then Vr 
should decrease in a large region of the disk, which 
should reduce the rate of matter actually accreted 
onto the black hole. And the average angular mo- 
mentum (< I >) of the disk increases with the 
increase of the peak and the width of the angular 
momentum distribution of the disk, which corre- 
sponds to the dips in emission, mass-loss param- 
eter and angular momentum loss rate. Although 
< I > oscillates with the same period as that of 
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shock, but the disk interestingly prefers to stay in 
the state < / > is greater than Imj- Since the 
disk itself is oscillating, all these flow parameters 
should oscillate with the same period. And indeed 
the bremsstrahlung emission, the mass loss rate, 
the angular momentum loss rate etc all oscillate 
with the same period of shock oscillation. 

4-2.2. Shock oscillation for higher viscosity 

The dynamics of the disk with higher viscosity 
parameter is different from that due to the lower 
one. For higher viscosity parameter the difference 
in the disk dynamics will arise from more efficient 
angular momentum transfer as well as higher vis- 
cous dissipation of heat, even if the outer bound- 
ary condition remains the same. In Figure 10, we 
have plotted the shock location with time (top) 
and the bremsstrahlung emission with time (bot- 
tom) for a fluid with the same injection param- 
eters as the inviscid flow described in §4.1, i.e., 
rj„j = lOOOrg, Vr{inj) = 2.970 x lO'^c, and 
Cs{inj) — 4.827 x lO^'^c and the viscosity pa- 
rameter is a — 0.1. The time variation of shock 
for a = 0.1 (Fig. 10) is distinctly different from 
that of a = 0.01 {i.e., Fig. 7). The inner shock 
forms, expands, and at some epoch collides with 
the contracting outer shock, while at some other 
epoch it disappears before colliding with the outer 
shock. The inner shock is weaker compared to the 
disk with lower a. The time evolution of the two 
shocks are somewhat similar to the initial phases 
of the shock variation for a = 0.01. Compari- 
son of the time variation of bremsstrahlung emis- 
sion with the time variation of the shock shows 
no correlation between shock minima and emis- 
sion maxima unlike the case for a = 0.01. In 
Figure 11a - lib, Vr (dashed-dot) and Cg (solid) 
are plotted corresponding to the emission max- 
ima (Fig. 11a) and emission minima (Fig. lib). 
Similarly the corresponding specific angular mo- 
mentum distributions are plotted for the emission 
maxima (Fig. 11c) and minima (Fig. lid), and 
the density too are plotted for the emission max- 
ima (Fig. lie) and minima (Fig. llf). The max- 
ima of the bremsstrahlung emission occurs when 
the inner shock is tending to form, while the min- 
ima occurs when the inner shock has not been 
formed (also refer Fig. 10). This change in the 
behavior of the shock and the emission properties 
compared to that of the a = 0.01, actually de- 



pends on the different rate of angular momentum 
transfer. Since the viscosity in the present case is 
ten fold higher than a = 0.01, the outward an- 
gular momentum transport is very efficient. So 
close to the black hole, the angular momentum 
rises steeply outward unlike the flow with lower 
viscosity {e.g., Fig. 11c - lid may be compared 
with Fig. 6a - 6d). If the shock is closer to the 
black hole then the peak of the angular momen- 
tum distribution is very close to the outer shock 
(Fig. lid); this causes matter to accrete more 
freely between the horizon and the peak of the an- 
gular momentum distribution, and hence the den- 
sity is lower (Fig. lib, lid and llf). This causes 
the emission to dip. As the shock moves out the 
angular momentum peak is situated further inside 
(Fig. 11c). And this causes the matter to madly 
fall inwards between the outer shock and the in- 
ner I peak. As the in falling matter encounters the 
angular momentum pile it decelerates drastically 
increasing the density considerably, and hence en- 
hances the bremsstrahlung emission (Fig. 11a, 
11c and lie). Eventually it forms an inner shock 
but the enhanced energy deposition in the post- 
inner shock region causes the inner shock to ex- 
pand and thereby reducing density and emission. 
In this connection one may point out that the 
immediate post-shock (for both inner and outer 
shocks) region may be decelerating or accelerating 
{e.g.. Fig s. 5a-5d, lla - llb). However, it was pre- 
dicted bv iNakavama ( 1992f) : Nobuta fc Hanawa 



( 19941 ) that post-shock acceleration and decelera- 
tion correspond to unstable and stable shocks, re- 
spectively. The reason for this is that no standing 
shock can exist in the viscous flow for the corre- 
sponding initial conditions. 

In the top panel of Figure 12, we plot M/Minj, 
and like the lower viscosity case its peak and 
trough coincides with that of the bremsstrahlung 
emission. The angular momentum loss rate 
L/Linj (middle) also follows the pattern of 
M/Minj. Since the angular momentum distri- 
bution is higher during the peak emission, the 
average angular momentum of the disk < I > 
(bottom) peak coincides with the emission peak. 
Moreover, the M/Minj < 1 most of the time, 
which means because of higher angular momen- 
tum most of the matter that are being injected in 
to the disk, are not being consumed by the black 
hole, which is indicated by the fact that < ^ > is 
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significantly higher than linj. 

The role of viscosity can be ascertained if one 
compares the viscous time with the advection time 
scale. The viscous time-scale may be defined as 

h J, 

-dr, (26) 

V 

where, v — p^j p and the advection time scale 

Tad = / (27) 

The time- variation of Tyis closely follows the cor- 
responding variation of the shock location, with 
minima and maxima of each, coinciding with the 
other at exactly the same time. When compared 



for a 



0.01 case, at the shock minima, Tv 



and Tad are comparable and we see the shock 
front is pushed outward. At the shock maxima 
Tad << TVis, i.e., advection dominates, conse- 
quently the shock front hurls inwards contract- 
ing significantly from few xlOO Vg to few xlO Vg. 
Hence the interplay between these two timescales, 
sustains the oscillation. The general pattern of 
the temporal behavior of the time scales and the 
relation to the shock oscillation for the flow with 
a = 0.1 is similar to that with the flow a = 0.01. 
However, Tad and Tvis for a — 0.1 are roughly 
comp arable and hence o s cillatio ns do not satu- 
rate. Foglizzo fc Tagger ( 2000f ) investigated for 
Bondi flow that entropy-acoustic cycles may sus- 
tain shock oscillation if Cs(rci)/cs(rs/i) ^ 1. In the 
simulations we have run, max \cs(rc.i)l Cs[rs-h)\ ~ 
10 but generally the ratio is less in most of the 
times. Hence the effect of entropy-acoustic cycles 
in regulating shock oscillation is probably moder- 
ate in our case. A multi-dimensional simulation 
may give a more definitive answer. 

5. SUMMARY AND DISCUSSION 

This paper is intended to study the time- 
dependent simulations of large amplitude oscil- 
lations of advective, viscous, sub-Keplerian disks, 
to complement earlier works of studying low am- 
plitude oscilla tions undertaken by Mo lteni and his 
collaborators ( Lanzafame et aZ]|l998f) . As an im- 
provement we have employed a new code which 
uses the Lagrangian TVD/remap approach. This 
code strictly conserved the angular momentum 



without viscosity, and reduced the numerical dis- 
sipation considerably (e.g., §3). Tests showed that 
the shock capturing capability of this code is bet- 
ter than both standard Eulerian code and La- 
grangian SPH code (e.g., Fig. 1), and followed 
the angular momentum transfer of the viscous, 
subsonic analytical solution extremely well (e.g., 
Fig. 2). 

Oscillation of accretion shock was borne out by 
the different rates of angular momentum transfer 
across the shock and the heat dissipated due to the 
presence of viscosity. It has been shown that in 
presence of low viscosity parameter (a — 0.003), 
the shock front of a disk, with the same initial 
and boundary conditions as those of the inviscid 
case, did tend to expand and settled at a larger 
distance from the disk (Fig. 4b). For an even 
higher viscosity (a ^ 0.005), the rate of angular 
momentum transfer was higher, which caused a 
faster rate of shock front expansion. As the shock 
front exceeded a possible equilibrium position it 
started to oscillate (Fig. 4c). However, it is 
to be remembered that the value of the critical 
viscosity parameter (a ~ 0.005 in the present 
case) is not sacrosanct, but actually depends on 
the initial condition. For example, it has been 
shown that the critical viscosity parameter will 
be higher for fiows with lower angular momen- 
tum, while for a fiuid with higher initial energy 
the critical viscosity para meter will be lower (see. 



Chakrabarti fc Das Il2004f) . Hence, if proper initial 



condition is used then a stable shock is expected to 
form for higher viscosity parameters (i.e., a ^ 0.1 
- 0.2) too, investigation of which, however, is not 
the point of interest for the present paper. 

A detailed study of the disk dynamics was con- 
ducted for reasonably high viscosity (i.e., a = 
0.01, & 0.1). For a = 0.01 the shock oscillation 
amplitude was found to be quite high ^ lOOr^. 
This resulted in a large sound speed gradient in 
the post-shock subsonic fiow. In case of large am- 
plitude shock oscillation, the rate of outward an- 
gular momentum transport in a region closer to 
the inner sonic point was shown to be much higher 
compared to the rate of angular momentum trans- 
port near the shock. As a result, our simulation 
showed that the angular momentum to be piled 
up in an intermediate region between the shock 
and the inner sonic point. The expanding shock 
also increased the inflow velocity in the immedi- 
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ate post-shock region only to be decelerated by the 
extra centrifugal pressure due to the piled-up an- 
gular momentum further inside the disk (e.g., Fig. 
6a - 6d). The inflow velocity in the post shock disk 
may be increased to the extent that it may again 
become supersonic, then the resistance from the 
excess centrifugal pressure from the piled-up angu- 
lar momentum distribution may cause the forma- 
tion of inner shock. In case of moderately high a, 
the distance between the peak of the angular mo- 
mentum distribution and the outer shock is large 
enough to allow for the Vr to become supersonic 
again and enhanced the possibility to form the in- 
ner shock. It is to be noted that the amplitude of 
shock oscillation will possibly be lesser for multi- 
dimensional simulation. Viscosity is more active in 
the post-shock disk, and hence the extra centrifual 
force due to the piled angular momentum and the 
heat dissipated by viscosity both actively take part 
in shock oscillation. However, in case of realistic 
accretion flow, a part of viscous heat dissipated in 
the post-shock disk will also be spend to puff it up, 
which would imply less outward push on the shock 
surface. Hence for a flow with same injection and 
viscosity parameters, the oscillation amplitude for 
a multi-dimensional disk is expected to be lesser 
compared to a purely conical flow. Consequently, 
the critical viscosity above which the disk becomes 
oscillatory will also be higher. 

The time evolution of shocks for higher viscos- 
ity was shown to be distinctly different from that 
of the lower one. The inner shock was weaker 
and more sporadic for a disk with a — 0.1. The 
main reason was because of the higher rate of an- 
gular momentum transport. Even when the shock 
was around lOOrg, highly efficient angular momen- 
tum transport created a smooth increase of angu- 
lar momentum, which only peaked closer to Vsh- 
As shock expanded Vr increased, but the opportu- 
nity to become supersonic was minimized since the 
peak of the /(r) was closer to the shock. Hence the 
inner shock, if formed at all, was weaker. However, 
since shock amplitude for a = 0.1 was much larger 
than the case with a — 0.01, with time the forma- 
tion of inner shock became more regular, and the 
behavior was more similar to that of a = 0.01. 

The oscillatory motion of the shock induced 
oscillation in all the disk parameters like emis- 
sion, rate of matter consumed by the black hole, 
the rate of angular momentum consumed by the 



black hole, and the average angular momentum 
of disk. All these parameters did oscillate with 
the same period as that of the shock. The disc 
oscillation started with a }t 0.005. Considering 
Mbh = 10 Mq, for a = 0.005 the oscillation 
frequency of the outer shock was 5 Hz, and the 
inner shock 10 Hz, for a — 0.006 the frequen- 
cies were 1 Hz and 3 Hz respectively, and for 
a = 0.01 the two frequencies were 0.125 Hz and 
0.25 Hz. Hence one may conclude that apart 
from the dependence of the oscillation frequency 
on injection parameters, the QPO frequency defi- 
nitely decreases with increasing viscosity and vice- 
versa. Observationally, GRO J1655-40 exhibits a 
rise in QPO frequency in its rising state and a 
fall i n QPO frequency in its declining phase in 
2005 (IChakrabarti et a^.boOSl) . IChakrabarti et all 
(|2009l ) plotted the QPO frequency for the object 



XTE J1550-564 in 1998 burst phase. They showed 
that in the rising phase of the outburst, the low 
frequency QPO increases from 0.08Hz to 13.1 Hz 
and then starts to decrease in the declining phase 
before disappearing. Such rise and fall of QPO fre- 
quencies may be explained by the change in shock 
oscillation frequency due to the change of the net 
viscosity of the disk. 

In presence of viscosity a positive angular mo- 
mentum gradient i.e., dl/dr > helps in outward 
transport of angular momentum. However, a neg- 
ative gradient may trigger inward transport of an- 
gular momentum. The dl/dr < condition was 
attained in the disk in at least two locations, at 
the outer shock front and just behind the peak 
of the specific angular momentum distribution. 
Those regions were subject to the rotational insta- 
bility, dl/dr < caused the average angular mo- 
mentum < I > oi the disk to increase, and hence 
the period and the amplitude of the shock oscil- 
lation to increase too. This is less perceptible for 
lower a and the shock oscillation achieved quasi- 
saturation, but for a = 0.1 the shock went outside 
the computation domain. We repeated the simu- 
lation with a = 0.3 (not presented in the paper) 
and in this case too the shock went outside the do- 
main, although formation of inner sonic point and 
oscillation of the two shocks were observed too. 

In case of multi-dimensional simulations, a part 
of the post shock matter would have ejected along 
the vertical direction in the form of winds, which 
would have carried away a part of the angular mo- 
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mentuni, such that the increase of < ^ > may have 
been arrested for higher a. This would have meant 
that the shock oscillation may saturate for a > 
0.1. Hence, wc conjecture that the non-saturation 
of shock oscillation for a ^ 0.1, could be an arti- 
fact of one-dimensional simulation. We will test it 
in a future work using multi-dimensional simula- 
tions. 
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Fig. 1. — Comparison of numerical solution with 
analytic solution for a onc-dinicinsional accretion 
flow in cylindrical geometry that allows a stand- 
ing shock. The plots show the density (p) and the 
radial Mach number (Mr). The solid curves rep- 
resent the analytical solution and the open circles 
represent the numerical solution with 2048 uni- 
form grid cells. 



This 2-column preprint was prepared with the AAS lAlJ^X 
macros v5.2. 
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Fig. 2.— Test of ADAF with 7 = 4/3 and a = 
0.3 under the Newtonian potential The sohd hnes 
represent the analytical self-similar solution, while 
the dots represent the numerical solution. The 
density p (a), and radial velocity ?v (lower curve 
of (b)) & adiabatic sound speed Cg (upper curve 
of (b)), specific angular momentum I (c), pressure 
p (d), are shown clockwise. 
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Fig. 3. — Comparison the analytical (solid) with 
the numerical solution (open circles) of hydrody- 
namical accretion shock for transonic flow with 
7 = 1.4 and I = l.Scrg. The injection ra- 
dius Tinj — lOOOrg, and at Vinj the radial ve- 
locity Vr{inj) = 2.970 x 10~^c and sound speed 
Csiinj) = 4.827 x lO'^c. 
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Fig. 4. — Comparison of the shock location rsh 
for a = 0.0 (a), a = 0.003 (b) and a = 0.006 
(c). The injection radius is rinj — lOOOr^, the 
injected radial velocity Vr{inj) = 2.970 x 10~^c, 
sound speed Cs{inj) = 4.827 x 10~^c, and angular 
momentum I 
7 = 1.4. 
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Fig. 5. — Four snapshots of radial velocity 
(dashed-dot) and sound speed (solid) of a vis- 
cous fluid of a = 0.01 and 7 = 1.4. The injec- 
tion parameters are rinj = lOOOrg, Imj = l.Scrg, 
Vr{inj) = 2.970 x lO-^c, and Cs{inj) = 4.827 x 
10^'^c. The time sequence goes as (a) — >■ (b) — >■ 
(c) ^ (d). 



(a) 1=2.900 xlOV{r/c) 


(b) t=-3.185 xlOV(ryc) 

J 


" 1 ' '1 

(c) t=3.400 xl0V{rg/c) 


1 1 

(d) t=3.615 xl0V(rg/G) 



10 100 1000 10 100 1000 



Fig. 6. — Specific angular momentum at the four 
different snapshots in the same simulation as in 
Fig. 5. The time sequence goes as (a) — >■ (b) — >■ 
(c) ^ (d). 
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Fig. 7. — Shock position (top) and bremsstrahlung 
emission (bottom) as a function of time in the 
same simulation as in Fig. 5. The upper curve 
in the top panel is the outer shock and the lower 
curve is the inner shock. The snap shots of time 
in Fig. 5a - 5d are marked on the top panel as a 
-d. 
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Fig. 8. — Left panels: Outer shock (top), 
inner shock (middle) and the bremsstrahlung 
emission (bottom) SiS Oil function on of time. 
Right panels: The power density spectra of the 
outer shock (top), inner shock (middle) and the 
bremsstrahlung emission (bottom). Same simula- 
tion as in Fig. 5. 
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Fig. 9. — Mass accretion rate (top), angular mo- 
mentum loss rate (middle) and averaged specific 
angular momentum (bottom) in the post shock re- 
gion of the outer shock as a function of time in the 
same simulation as in Fig. 5. 
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Fig. 10. — Shock position (top) and 

bremsstrahlung emission (bottom) as a func- 
tion of time. The upper curve in the top panel 
is the outer shock and the lower curve is the 
inner shock. The simulation is for a = 0.1 and 
7 = 1.4, and injection parameters rj„j = lOOOrg, 
linj = l.Scrg, Vr{inj) = 2.970 x 10~^c, and 
Cs,{inj) = 4.827 x lO'^c. 



Fig. 11. — Radial velocity (dashed-dot) and 
sound speed (solid) of viscous fluid for time t = 
2.65x10^3 (a) and t = 2.87xlO-Vg (b). Spe- 
cific angular momentum is plotted for time t = 
2.65xlOVg (c) and t = 2.87x10^X9 (d). Density 
in arbitrary units is plotted for t = 2.65xlO^Tg 
(e) and t = 2.87xlOVg (f). The simulation is 
for a = 0.1, and 7 = 1.4. The injection param- 
eters are Vinj = lOOOr^, linj = IScrg, Vr{inj) = 
2.970 x lO-^c, and Cs{inj) = 4.827 x lO'^c. 
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Fig. 12. — Mass accretion rate (top), angular mo- 
mentum loss rate (middle) and averaged specific 
angular momentum (bottom) in the post shock re- 
gion of the outer shock as a function of time in the 
same simulation as in Fig. 10. 
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